## -------------------------------------------------------------------------- ##
##  "Who would defend Kinmen and Matsu (and what is Wuchiu?):
##   Taiwanese public opinion on outlying islands"
##  Authors: Andrew Chubb, Ronan Tse-min Fu, Elaine I-lien Lee
##
##  This file reproduces Figures 3-12 in the paper.
## -------------------------------------------------------------------------- ##

sessionInfo() # R version 4.3.3

rm(list = ls())
setwd("~/R/CFL_replication_files/")


## ---- Install and load required packages -------------------------------------

pkgs <- c("Rmisc", "tidyverse", "estimatr", "marginaleffects", "patchwork")
for (pkg in pkgs) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg)
  }
  library(pkg, character.only = TRUE)
}


## ---- Data loading & recoding ------------------------------------------------

## 2023 survey (N = 1,541)
load("df2023.rda")

## 2024 survey + experiment (N = 1,403)
load("df2024.rda")


## ---- Derived objects --------------------------------------------------------

## df2023_long: long format for mean and OLS plots (Fig 3)
df2023_long <- df2023 %>%
  rowid_to_column("ID") %>% 
  pivot_longer(cols = c(Taiping, Kinmen, Penghu, Matsu),
               names_to  = "island",
               values_to = "defending") %>%
  mutate(island = factor(island, levels = c("Penghu", "Kinmen", "Matsu", "Taiping")))


## df2024_long: long format for belonging/security battery (Fig 4, 10, 12)
island_order <- c("Penghu", "Kinmen", "Matsu", "Taiping", "Wuchiu", "Pratas")

## Place name mapping: suffix 1=Kinmen, 2=Matsu, 3=Wuchiu, 4=Taiping, 5=Pratas, 6=Penghu
place_mapping <- tibble(
  variable_suffix = 1:6,
  place_name = c("Kinmen", "Matsu", "Wuchiu", "Taiping", "Pratas", "Penghu")
)

df2024_long <- df2024 %>%
  rowid_to_column("ID") %>% 
  pivot_longer(
    cols      = starts_with("belong_") | starts_with("security_"),
    names_to  = "variable",
    values_to = "response"
  ) %>%
  mutate(
    question_type   = case_when(
      grepl("belong",   variable) ~ "Belong to Taiwan",
      grepl("security", variable) ~ "Important to Taiwan's security"
    ),
    variable_suffix = as.integer(substr(variable, nchar(variable), nchar(variable)))
  ) %>%
  left_join(place_mapping, by = "variable_suffix") %>%
  select(ID, variable, response, question_type, place_name, everything()) %>% 
  mutate(response_bi = case_when(response == "Yes" ~ 1,
                                 response == "No" ~ 0,
                                 TRUE ~ NA)) %>% 
  select(-variable_suffix, -starts_with("exp"))

rm(place_mapping)


## ---- Shared theme -----------------------------------------------------------

theme_all <- theme(
  panel.background   = element_blank(),
  axis.line          = element_line(color = "black"),
  axis.title.y       = element_text(size = 24, color = "black", margin = margin(r = 8, l = 6), face = "bold"),
  axis.text.y        = element_text(size = 17, color = "black", margin = margin(r = 3)),
  axis.text.x        = element_text(size = 22, color = "black", margin = margin(t = 5, b = 2)),
  axis.ticks.x       = element_blank(),
  strip.text         = element_text(size = 20, color = "black", face = "bold", margin = margin(t = 8, b = 8, r = 8, l = 8)),
  panel.grid.major.y = element_line(color = "gray70", linewidth = .25, linetype = "dotted"),
  legend.position    = "none"
)


## ========================================================================== ##
## FIGURE 3 -----
## 2023 survey: mean approval of defending four islands (0-10 scale)
## Data: df2023_long
## x: island (Penghu / Kinmen / Matsu / Taiping)
## ========================================================================== ##

df2023_long %>%
  group_by(island) %>%
  summarize(mean  = mean(defending, na.rm = TRUE),
            ci    = as.numeric(t.test(defending)$conf.int[2] - t.test(defending)$conf.int[1]) / 2,
            lower = mean - ci,
            upper = mean + ci) %>%
  ungroup() %>%
  ggplot(aes(y = mean, x = island, color = island)) +
  geom_hline(aes(yintercept = 5), linetype = "solid", color = "gray90") +
  geom_pointrange(aes(ymin = lower, ymax = upper), lwd = 1.9, size = 1.2) +
  scale_y_continuous(limits = c(4, 8), breaks = seq(4, 8, 1)) +
  scale_x_discrete(limits = c("Penghu", "Kinmen", "Matsu", "Taiping")) +
  labs(y = "Average Approval of Defending", x = NULL) +
  theme_all


## ========================================================================== ##
## FIGURE 4 -----
## 2024 survey: mean binary agreement — belong to Taiwan vs. important to security
## Data: df2024_long
## DV: response_bi (0/1); facet: question_type; x: place_name (6 islands)
## ========================================================================== ##

df2024_long %>%
  group_by(question_type, place_name) %>%
  summarize(mean  = mean(response_bi, na.rm = TRUE),
            ci    = as.numeric(t.test(response_bi)$conf.int[2] - t.test(response_bi)$conf.int[1]) / 2,
            lower = mean - ci,
            upper = mean + ci) %>%
  ungroup() %>%
  mutate(place_name = factor(place_name, levels = island_order)) %>%
  ggplot(aes(y = mean, x = place_name, color = place_name)) +
  facet_grid(. ~ question_type) +
  geom_pointrange(aes(ymin = lower, ymax = upper), lwd = 1.9, size = 1.2) +
  scale_y_continuous(limits = c(0.8, 1), breaks = seq(.8, 1, .05),
                     labels = scales::percent_format(suffix = "%")) +
  scale_x_discrete(limits = island_order) +
  labs(y = "Average Agreement that...", x = NULL) +
  theme_classic() +
  theme_all +
  theme(axis.text.x     = element_text(size = 20, color = "black", margin = margin(t = 5, b = 2)),
        panel.spacing.x = unit(.5, "cm"))


## ========================================================================== ##
## FIGURE 5 -----
## Experiment: ATE of territory (exp_t2)
## Data: df2024
## Contrasts: Wuchiu - Taiwan, Matsu - Taiwan  (base: Taiwan)
## Model: lm_robust(exp_defend ~ exp_t1 * exp_t2)
## ========================================================================== ##

tidy(avg_comparisons(lm_robust(exp_defend ~ exp_t1 * exp_t2, data = df2024))) %>%
  filter(term == "exp_t2") %>%
  ggplot(aes(x = estimate, y = contrast, color = contrast)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "darkgray") +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), lwd = 1.8, size = 1.2) +
  scale_x_continuous(limits = c(-.23, .23), breaks = seq(-.2, .2, .1)) +
  scale_y_discrete(limits = rev(c("Territory", "Matsu - Taiwan", "Wuchiu - Taiwan")),
                   labels = rev(c("Territory (base: Taiwan)", "        Matsu", "        Wuchiu"))) +
  labs(x = "Average Effect on Approval of Defending", y = NULL) +
  theme(panel.background = element_blank(),
        axis.line        = element_line(color = "black"),
        axis.title.x     = element_text(size = 24, color = "black", margin = margin(t = 8, b = 3), face = "bold"),
        axis.text.x      = element_text(size = 17, color = "black", margin = margin(t = 3)),
        axis.text.y      = element_text(size = 21, color = "black", margin = margin(r = 2, l = 2), hjust = 0),
        axis.ticks.y     = element_blank(),
        legend.position  = "none")


## ========================================================================== ##
## FIGURE 6 -----
## Experiment: ATE of attack type (exp_t1)
## Data: df2024
## Contrast: Seizure - Blockade  (base: Blockade)
## Model: lm_robust(exp_defend ~ exp_t1 * exp_t2)
## ========================================================================== ##

tidy(avg_comparisons(lm_robust(exp_defend ~ exp_t1 * exp_t2, data = df2024))) %>%
  filter(term == "exp_t1") %>%
  ggplot(aes(x = estimate, y = contrast)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "darkgray") +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), lwd = 1.8, size = 1.2) +
  scale_x_continuous(limits = c(-.23, .23), breaks = seq(-.2, .2, .1)) +
  scale_y_discrete(limits = rev(c("Chinese attack", "Seizure - Blockade")),
                   labels = rev(c("Chinese attack (base: Blockade)", "        Seizure"))) +
  labs(x = "Average Effect on Approval of Defending", y = NULL) +
  theme(panel.background = element_blank(),
        axis.line        = element_line(color = "black"),
        axis.title.x     = element_text(size = 24, color = "black", margin = margin(t = 8, b = 3), face = "bold"),
        axis.text.x      = element_text(size = 17, color = "black", margin = margin(t = 3)),
        axis.text.y      = element_text(size = 21, color = "black", margin = margin(r = 2, l = 2), hjust = 0),
        axis.ticks.y     = element_blank(),
        legend.position  = "none")


## ========================================================================== ##
## FIGURE 7 -----
## Experiment: predicted means by 2x3 treatment cell
## Data: df2024
## DV: exp_defend (1-5); facet: exp_t1; x: exp_t2
## Model: lm_robust(exp_defend ~ exp_t1 * exp_t2)
## ========================================================================== ##

tidy(avg_predictions(
  lm_robust(exp_defend ~ exp_t1 * exp_t2, data = df2024),
  by = c("exp_t1", "exp_t2")
  )) %>%
  ggplot(aes(y = estimate, x = exp_t2, color = exp_t2)) +
  facet_grid(. ~ exp_t1) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), lwd = 1.9, size = 1.2) +
  scale_y_continuous(limits = c(3, 4), breaks = seq(3, 4, .2)) +
  scale_x_discrete(limits = c("Taiwan", "Wuchiu", "Matsu")) +
  labs(y = "Average Approval of Defending", x = NULL) +
  theme_classic() +
  theme_all


## ========================================================================== ##
## FIGURE 8 -----
## 2023 survey: subgroup means — 6-panel
## Data: df2023
## DV: defending (0-10); x: island (Taiping / Kinmen / Penghu / Matsu)
## Panels: gender / age / edu / identity / party / tondu
## ========================================================================== ##

{
  measure_vars <- c("Taiping", "Kinmen", "Penghu", "Matsu")
  group_vars   <- c("gender", "age", "edu", "identity", "party", "tondu")
  all_results  <- list()
  
  for (group in group_vars) {
    results <- data.frame()
    for (measure in measure_vars) {
      temp        <- Rmisc::summarySE(df2023, measurevar = measure, groupvar = group, na.rm = TRUE)
      temp$island <- measure
      results     <- bind_rows(results, temp)
      rm(temp)
    }
    results <- results %>%
      mutate(fight_if = coalesce(Taiping, Kinmen, Penghu, Matsu)) %>%
      select(-c(Taiping, Kinmen, Penghu, Matsu))
    all_results[[group]] <- results
    rm(results)
  }
  
  dat_gender   <- all_results[["gender"]]
  dat_age      <- all_results[["age"]]
  dat_edu      <- all_results[["edu"]]
  dat_identity <- all_results[["identity"]]
  dat_party    <- all_results[["party"]]
  dat_tondu    <- all_results[["tondu"]]
  rm(all_results, measure_vars, group_vars)
}

plot_layer_sub <- list(
  geom_hline(aes(yintercept = 5), linetype = "solid", color = "gray90"),
  geom_point(aes(shape = sub), size = 3.7, position = position_dodge(.4)),
  geom_linerange(aes(ymin = lower, ymax = upper), lwd = 1.4,
                 position = position_dodge(.4), show.legend = FALSE),
  scale_y_continuous(limits = c(3.25, 8.75), breaks = seq(4, 8, 1)),
  scale_shape_manual(values = c(15, 17, 19, 18)),
  labs(y = "Average Approval of Defending", x = NULL),
  theme(panel.background   = element_blank(),
        plot.title         = element_text(size = 22, color = "black",
                                          margin = margin(t = 10, b = 8), face = "bold", hjust = 0.5),
        axis.line          = element_line(color = "black"),
        axis.title.y       = element_text(size = 19, color = "black", margin = margin(r = 10), face = "bold"),
        axis.text.y        = element_text(size = 14, color = "black", margin = margin(r = 3)),
        axis.text.x        = element_text(size = 16, color = "black", margin = margin(t = 4, b = 4)),
        axis.ticks.x       = element_blank(),
        panel.grid.major.y = element_line(color = "gray70", linewidth = .25, linetype = "dotted"),
        legend.text        = element_text(size = 16),
        legend.title       = element_blank(),
        legend.key.spacing.x = unit(.3, "cm"),
        legend.position    = "top")
)

dat_gender %>% dplyr::rename(sub = gender) %>% filter(sub != "") %>%
  mutate(lower = fight_if - ci, upper = fight_if + ci) %>%
  ggplot(aes(y = fight_if, x = island, color = sub)) + plot_layer_sub + labs(title = "By Gender") +

  dat_age %>% dplyr::rename(sub = age) %>%
  mutate(lower = fight_if - ci, upper = fight_if + ci) %>%
  ggplot(aes(y = fight_if, x = island, color = sub)) + plot_layer_sub + labs(title = "By Age") +
  
  dat_edu %>% dplyr::rename(sub = edu) %>%
  mutate(lower = fight_if - ci, upper = fight_if + ci) %>%
  ggplot(aes(y = fight_if, x = island, color = sub)) + plot_layer_sub + labs(title = "By Education") +
  
  dat_identity %>% dplyr::rename(sub = identity) %>% filter(sub != "") %>%
  mutate(lower = fight_if - ci, upper = fight_if + ci) %>%
  ggplot(aes(y = fight_if, x = island, color = sub)) + plot_layer_sub + labs(title = "By Ethnic ID") +
  
  dat_party %>% dplyr::rename(sub = party) %>% filter(sub != "") %>%
  mutate(lower = fight_if - ci, upper = fight_if + ci) %>%
  mutate(sub = factor(sub, levels = c("KMT", "TPP", "DPP"))) %>%
  ggplot(aes(y = fight_if, x = island, color = sub)) + plot_layer_sub + labs(title = "By Party ID") +
  
  dat_tondu %>% dplyr::rename(sub = tondu) %>% filter(sub != "") %>%
  mutate(lower = fight_if - ci, upper = fight_if + ci) %>%
  mutate(sub = factor(sub, levels = c("Unification", "Status Quo", "Independence"))) %>%
  ggplot(aes(y = fight_if, x = island, color = sub)) + plot_layer_sub + labs(title = "By Tondu") +
  
  plot_layout(ncol = 3, axis_titles = "collect")


## ========================================================================== ##
## FIGURE 9 -----
## Experiment: predicted approval by identity (Taiwanese vs. Dual/Chinese)
## Data: df2024, subsetted by identity
## Facet: exp_t1 (Blockade / Seizure); x: exp_t2
## identity: Taiwanese / Dual/Chinese
## ========================================================================== ##

rbind(
  tidy(avg_predictions(
    lm_robust(exp_defend ~ exp_t1 * exp_t2, data = subset(df2024, identity == "Dual/Chinese")),
    by = c("exp_t1", "exp_t2"))) %>% mutate(identity = "Dual/Chinese"),
  tidy(avg_predictions(
    lm_robust(exp_defend ~ exp_t1 * exp_t2, data = subset(df2024, identity == "Taiwanese")),
    by = c("exp_t1", "exp_t2"))) %>% mutate(identity = "Taiwanese")
) %>%
  ggplot(aes(y = estimate, x = exp_t2, color = identity, shape = identity)) +
  facet_grid(. ~ exp_t1) +
  geom_point(size = 4.5, position = position_dodge(.5)) +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), lwd = 1.8,
                 position = position_dodge(.5), show.legend = FALSE) +
  scale_y_continuous(limits = c(2.4, 4.6), breaks = seq(2.5, 4.5, .5)) +
  scale_x_discrete(limits = c("Taiwan", "Wuchiu", "Matsu")) +
  labs(y = "Average Approval of Defending", x = NULL) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf) +
  theme_classic() +
  theme_all +
  theme(legend.text          = element_text(size = 19),
        legend.title         = element_blank(),
        legend.key.spacing.x = unit(.8, "cm"),
        legend.position      = "top")


## ========================================================================== ##
## FIGURE 10 -----
## 2024 survey: belong vs. security by identity (Taiwanese vs. Dual/Chinese)
## Data: df2024_long
## DV: response_bi (0/1); facet: question_type; x: place_name
## identity: Taiwanese / Dual/Chinese
## ========================================================================== ##

df2024_long %>%
  na.omit(identity) %>% 
  group_by(question_type, place_name, identity) %>%
  summarize(mean  = mean(response_bi, na.rm = TRUE),
            ci    = as.numeric(t.test(response_bi)$conf.int[2] - t.test(response_bi)$conf.int[1])/2,
            lower = mean - ci,
            upper = mean + ci) %>% 
  ungroup() %>% 
  mutate(place_name = factor(place_name, levels = island_order)) %>%
  ggplot(aes(y = mean, x = place_name, color = identity, shape = identity)) +
  facet_grid(. ~ question_type) +
  geom_point(size = 4.5, position = position_dodge(.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), lwd = 1.8,
                 position = position_dodge(.5), show.legend = FALSE) +
  scale_y_continuous(limits = c(0.74, 1.05), breaks = seq(.8, 1, .1)) +
  scale_x_discrete(limits = island_order) +
  labs(y = "Average Agreement that...", x = NULL) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf) +
  theme_classic() +
  theme_all +
  theme(axis.text.x          = element_text(size = 20, color = "black", margin = margin(t = 5, b = 2)),
        panel.spacing.x      = unit(.4, "cm"),
        legend.text          = element_text(size = 19),
        legend.title         = element_blank(),
        legend.key.spacing.x = unit(.8, "cm"),
        legend.position      = "top")


## ========================================================================== ##
## FIGURE 11 -----
## Experiment: predicted approval by identity2 (Chinese vs. Dual/Taiwanese)
## Data: df2024, subsetted by identity2
## Facet: exp_t1 (Blockade / Seizure); x: exp_t2
## identity2: Chinese / Dual/Taiwanese
## ========================================================================== ##

rbind(
  tidy(avg_predictions(
    lm_robust(exp_defend ~ exp_t1 * exp_t2, data = subset(df2024, identity2 == "Chinese")),
    by = c("exp_t1", "exp_t2"))) %>% mutate(identity = "Chinese"),
  tidy(avg_predictions(
    lm_robust(exp_defend ~ exp_t1 * exp_t2, data = subset(df2024, identity2 == "Dual/Taiwanese")),
    by = c("exp_t1", "exp_t2"))) %>% mutate(identity = "Dual/Taiwanese")
) %>%
  ggplot(aes(y = estimate, x = exp_t2, color = identity, shape = identity)) +
  facet_grid(. ~ exp_t1) +
  geom_point(size = 4.5, position = position_dodge(.5)) +
  geom_linerange(aes(ymin = conf.low, ymax = conf.high), lwd = 1.8,
                 position = position_dodge(.5), show.legend = FALSE) +
  scale_y_continuous(limits = c(0, 4.4), breaks = seq(0, 4, 1)) +
  scale_x_discrete(limits = c("Taiwan", "Wuchiu", "Matsu")) +
  labs(y = "Average Approval of Defending", x = NULL) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf) +
  theme_classic() +
  theme_all +
  theme(legend.text          = element_text(size = 19),
        legend.title         = element_blank(),
        legend.key.spacing.x = unit(.8, "cm"),
        legend.position      = "top")


## ========================================================================== ##
## FIGURE 12 -----
## 2024 survey: belong to Taiwan only, by identity2 (Chinese vs. Dual/Taiwanese)
## Data: df2024_long, filtered to question_type == "Belong to Taiwan"
## DV: response_bi (0/1); x: place_name
## identity2: Chinese / Dual/Taiwanese
## ========================================================================== ##

df2024_long %>%
  na.omit(identity2) %>% 
  group_by(question_type, place_name, identity2) %>%
  summarize(mean  = mean(response_bi, na.rm = TRUE),
            ci    = as.numeric(t.test(response_bi)$conf.int[2] - t.test(response_bi)$conf.int[1]) / 2,
            lower = mean - ci,
            upper = mean + ci) %>%
  ungroup() %>% 
  filter(question_type == "Belong to Taiwan") %>% 
  mutate(place_name = factor(place_name, levels = island_order)) %>%
  ggplot(aes(y = mean, x = place_name, color = identity2, shape = identity2)) +
  geom_point(size = 4.5, position = position_dodge(.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), lwd = 1.8,
                 position = position_dodge(.5), show.legend = FALSE) +
  scale_y_continuous(limits = c(0.58, 1.03), breaks = seq(.6, 1, .1)) +
  scale_x_discrete(limits = island_order) +
  labs(y = "Average Agreement that Belong to Taiwan", x = NULL) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf) +
  theme_classic() +
  theme_all +
  theme(axis.text.x          = element_text(size = 20, color = "black", margin = margin(t = 5, b = 2)),
        legend.text          = element_text(size = 19),
        legend.title         = element_blank(),
        legend.key.spacing.x = unit(.8, "cm"),
        legend.position      = "top")

